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We address the dynamical and statistical description of stably stratified turbulent boundary layers 
with the important example of the atmospheric boundary layer with a stable temperature stratifi- 
cation in mind. Traditional approaches to this problem, based on the profiles of mean quantities, 
velocity second-order correlations, and dimensional estimates of the turbulent thermal flux run into 
a well known difficulty, predicting the suppression of turbulence at a small critical value of the 
Richardson number, in contradiction with observations. Phenomenological attempts to overcome 
this problem suffer from various theoretical inconsistencies. Here we present an approach taking 
into full account all the second-order statistics, which allows us to respect the conservation of total 
mechanical energy. The analysis culminates in an analytic solution of the profiles of all mean quan- 
tities and all second-order correlations removing the unphysical predictions of previous theories. We 
propose that the approach taken here is sufficient to describe the lower parts of the atmospheric 
boundary layer, as long as the Richardson number does not exceed an order of unity. For much 
higher Richardson numbers the physics may change qualitatively, requiring careful consideration of 
the potential Kelvin-Helmoholtz waves and their interaction with the vortical turbulence. 
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A 


Thermal flux production vector, (1.7c) 





Buoyancy parameter, g(3 


B 


Pressure-temperature-gradient-vector, (1.7f) 




Thermal expansion coefficient, (A6d) 


c 


Energy conversion tensor, (1.7b) 


7m 


Relaxation frequency of T\j , i = j 


v/vt 


Substantial derivative, d/dt + W-V 


7ri 


Relaxation frequency of , i =/= j 


D/Dt 


Mean substantial derivative, d/dt + U ■ V 


7rd 


Relaxation frequency of F 


E K 


Turbulent kinetic energy per unit mass, \\u\ 2 


iu u 


Relaxation frequency of E K 


E e 


"Temperature energy" per unit mass, ^9 2 


lee 


Relaxation frequency of E e 


F 


Turbulent thermal flux per unit mass, (u9) 




Dissipation tensor of Ty, (1.6) 


F* 


Thermal flux at zero elevation z = 


e 


Dissipation vector of F, (1.6) 


g 


Gravity acceleration, g = —gz 


£ 


Dissipation of E e , (1.6) 


L 


Monin-Obukhov length, ul//3F* 


e 


Total potential temperature 


I 


Outer scale of turbulence, external parameter 


e d 


Deviation of from BRS 


rij 


Rate of Reynolds stress production, (1.7a) 


e 


Mean potential temperature, (@d) 


P, P, V* 


Total, fluctuating and zero level pressures 





Fluctuating potential temperature 


Pr T 


Turbulent Prandtl number, v T /Xi 


o,„ 


Potential temperature 


Riflux 


Flux Richardson number, (3F z /t xz S v 




at zero elevation, F*/u* 


Rigrad 


Gradient Richardson number, f3S /S 2 


A* 


Viscous lengthscale, u/u* 


Ttij 


Pressure-rate-of-strain-tensor, (1.7c) 


V 


Kinematic viscosity 


Sjj 


Mean velocity gradient, dU/dz 


;/ T 


Turbulent viscosity 


iS e 


Mean potential temperature gradient, dQ/dz 


P 


density of the fluid 


T 


Molecular temperature 


T u 


Reynolds stress tensor, (uiUj) 


U 


Velocity field 


T* 


Mechanical momentum flux 


u 


Mean velocity, (U) 




at zero elevation (at the graund) 


u 


Fluctuating velocity, U — U 


X 


dynamical thermal conductivity 


it* 


(Wall) friction velocity, ^frl 


Xt 


Turbulent thermal conductivity 


X, z 


stream-vise and vertical (wall-normal) directions 


BRS 


Basic Reference State 



Introduction spheric boundary layer, this is the part of the atmosphere 

where the surface influences the temperature, moisture, 

The lower levels of the atmosphere are usually strongly 
influenced by the Earth's surface. Known as the atmo- 
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and velocity of the air above through the turbulent trans- 
fer of mass. 

The stability of the atmospheric boundary layer de- 
pends on the profiles of the density and the temperature 
as a function of the height above the ground. During 
normal summer days the land mass warms up and the 
temperature is higher at lower elevations. If it were not 
for the decrease in density of the air as a function of 
the height, such a situation of heating from below would 
have been always highly unstable. In fact, the boundary 
layer is considered stable as long as the temperature de- 
creases at the dry adiabatic lapse rate (T" w — 9.8°C per 
kilometer) throughout most of the boundary layer. With 
such a rate of cooling one balances out the decrease in 
density. With a higher degree of cooling one refers to 
the atmospheric boundary layer as unstably stratified, 
whereas with a lower degree of cooling the situation is 
stably stratified. Stably stratified boundary layer oc- 
curs typically during clear, calm nights. In extreme cases 
turbulence tends to cease, and radiational cooling from 
the surface results in a temperature that increases with 
height above the surface. 

The tendency of the atmosphere to be turbulent does 
not depend only on the rate of cooling but also on the 
mean shear in the vertical direction. The commonly used 
parameter to describe the tendency of the atmosphere 
to be turbulent is the "gradient" Richardson number 
(Richardson, 1920), defined as 



the case of thermal stratification): 

nj = (uiUj) , F = (u9) 



(0.2) 



Ri. 



grad 



(3dQ(z)/dz 
[dU x /dz] 2 



(0.1) 



where x is the stream- wise direction, z is the height above 
the ground, Q(z) is the the "mean potential tempera- 
ture" profile, (which differs from the mean temperature 
profile T(z) by accounting for the adiabatic cooling of the 
air during its expansion: d@(z)/dz = dT{z)/dz + \T'\), 
[3 — (3g is the buoyancy parameter in which (3 is the adi- 
abatic thermal expansion coefficient, [defined below in 
Eq. (A6d); for an ideal gas j3 = 1/T], and g is the gravi- 
tational acceleration. The mean shear dU x /dz is defined 
in terms of the mean velocity U, which in the simplest 
case of flat geometry depends only on the vertical coordi- 
nate z. The parameter Ri gra d represents the ratio of the 
generation or suppression of turbulence by buoyant pro- 
duction of energy to the mechanical generation of energy 
by wind shear. 

In this paper we consider the description of stably 
stratified turbulent boundary layers (TBL), taking as an 
example the case of stable thermal stratification. Since 
the 50's of twentieth century, traditional models of strat- 
ified TBL generalize models of unstratified TBL, based 
on the budget equations for the kinetic energy and me- 
chanical momentum; see reviews of Umlauf and Burchard 
(2005), Weng and Taylor (2003). The main difficulty is 
that the budget equations are not closed; they involve 
turbulent fluxes of mechanical moments (known as 
the "Reynolds stress" tensor) and a thermal flux F (for 



where u and 9 stand for the turbulent fluctuating velocity 
and the potential temperature with zero mean. The na- 
ture of the averaging procedure behind the symbol (• • • ) 
will be specified below. 

Earlier estimates of the fluxes (0.2) are based on the 
concept of the down-gradient turbulent transport, in 
which, similarly to the case of molecular transport, the 
flux is taken proportional to the gradient of transported 
property times a corresponding (turbulent) transport co- 
efficient: 

t xz = -v T dU x /dz, v T w C v l zS fr z ~ z , (0.3a) 
F z = -XrdQ/dz, Xt-C x £ zV ^T z , etc. (0.3b) 

Here the turbulent-eddy viscosity i/ T and turbulent ther- 
mal conductivity % T are estimated by dimensional rea- 
soning via the vertical turbulent velocity y/r^ and a scale 
l z (which in the simplest case is determined by the ele- 
vation z). The dimcnsionless coefficients C v and C x are 
assumed to be of the order of unity. 

This approach meets serious difficulties (Zcman, 1981), 
in particular, it predicts a full suppression of turbulence 
when the stratification exceeds a critical level, for which 
Rigrad =Ri C r ~ 0.25. On the other hand, in observa- 
tions of the atmospheric turbulent boundary layer tur- 
bulence exists for much larger values than Ri gra d = 0.25: 
experimentally above Ri gra d = 10 and even more (for 
detailed discussion see Zilitinkevich et al, submitted to 
Science). In models for weather predictions this prob- 
lem is "fixed" by introducing fit functions C„(Ri gra d) 
and C x (Rigrad) instead of the constant C v and C x in the 
model parametrization (0.3). This technical "solution" 
is not based on a physical derivation and just masks the 
shortcomings of the model. To really solve the problem 
one has to understand its physical origin, even though 
from a purely formal viewpoint it is indeed possible that 
a dimensionless coefficient like C x can be any function of 

Rigrad • 

To expose the physical reason for the failure of the 
down-gradient approach, recall that in a stratified flow, 
in the presence of gravity, the turbulent kinetic energy 
is not an integral of motion. Only the total mechanical 
energy, the sum of the kinetic and the potential energy, is 
conserved in the inviscid limit. As it was shown already 
by Richardson, the difficult point is that an important 
contribution to the potential energy comes not just from 
the mean density profile, but from the density fluctu- 
ations. Clearly, any reasonable model of the turbulent 
boundary layer must obey the conservations laws. 

The physical requirement of conserving the total me- 
chanical energy calls for an explicit consideration not 
only of the mean profiles, but also of all the relevant 
second-order, one-point, simultaneous correlation func- 
tions of all the fluctuating fields together with some clo- 
sure procedure for the appearing third order moments. 
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First of all, in order to account for the important ef- 
fect of stratification on the anisotropy, we must write 
explicit equations for the entire Reynolds stress tensor, 
Tij = (uiUj) . Next, in the case of the temperature strat- 
ified turbulent boundary layers we follow tradition [see, 
e.g. see e.g. Zeeman, (1981), Hunt et al. (1988), Schu- 
mann and Gerz (1995), Hanazaki and Hunt (2004), Keller 
and van Atta (2000), Stretch et al. (2001), Elperin et al. 
(2002), Cheng et al. (2002) Luytcn et al. (2002), and 
Rchmann and Hwang (2005)] and account for the turbu- 
lent potential energy which is proportional to the vari- 
ance of the potential temperature deviation, (9 2 ). And 
last but not least, we have to consider explicitly equa- 
tions for the the vertical fluxes, t xz and F z , which include 
the down gradient terms proportional to the velocity and 
temperature gradients, and counter-gradient terms, pro- 
portional to F x (in the equation for t xz ) and to (0 2 ) (in 
the equation for F z ) . 

Unfortunately, the resulting second order closure seems 
to be inconsistent with the variety of boundary-layer 
data, and many authors took the liberty to introduce 
additional fitting parameters and sometimes fitting func- 
tions to achieve a better agreement with the data (see 
reviews of Umlauf and Burchard (2005), Wcng and Tay- 
lor (2003), Zeeman, (1981), Melor and Yamada (1974), 
and references therein). Moreover, in the second order 
closures the problem of critical Richardson number, Ri cr , 
seems to persists (Cheng et al., 2002; Canuto, 2002). 

An interesting attempt to improve the modeling of sta- 
bly stratified flows with the aim to remove Ri cr is re- 
ported in a paper by Zilitinkevich, Elperin, Kleeorin and 
Rogachevskii (2007) in this issue, referred hereafter to as 
the ZEKR-paper. Following the traditional second-order 
scheme (see e.g. Zeeman, 1981), they removed the un- 
warranted approximation (0.3b), and replace it with an 
exact equation for the vertical thermal flux F z and for 
the "temperature energy" (6 2 ^ /2, using dimensional rea- 
soning to estimate required third-order correlations. To 
account for the important effect of stratification on the 
anisotropy of turbulence, they, in agreement with Cheng 
et al. (2002), explicitly considered budget equations for 
the partial kinetic energies Tu/2. Nevertheless ZEKR 
do not close rigorously the momentum budget equation. 
They combine three terms in this equations into, what 
they called "effective dissipation rate" , arriving actually 
again at the down-gradient approximation (0.3a) for the 
vertical flux of x-componcnt of the mechanical moment, 
t xz . In our paper we treat these three terms separately, 
considered also explicitly the horizontal heat flux F x . 

Notice that in spite of obvious inconsistency of the 
first-order schemes, "most of the practically used tur- 
bulent models are based on the concept of the down- 
gradient transport". One of the reasons is that in the 
second-order schemes instead of two down-gradient equa- 
tions (0.3) one needs to take into account eight nonlin- 
ear coupled additional equations i.e. four equations for 
the Reynolds stresses, three equations for the heat fluxes 
and equation for the temperature variance. As the result, 



second-order schemes have seemed to be rather cumber- 
some for comprehensive analytical treatment and have 
allowed to find only some relationships between correla- 
tion functions (see, e.g., Cheng et al., 2002). Unfortu- 
nately, the numerical solutions to the complete set of the 
second-order schemes equations which involve too many 
fitting parameters are much less informative in clarifica- 
tion of physical picture of the phenomenon than desired 
analytical ones. 

In this paper we suggest a relatively simple second- 
order closure model of turbulent boundary layer with sta- 
ble temperature stratification that, from the one hand, 
accounts (as we believe) for main relevant physics in 
the stratified TBL and, from the other hand, is simple 
enough to allow complete analytical treatment including 
the problem of critical Ri gra d- To reach this goal we ap- 
proximate the third order correlations via the first- and 
second-order ones, accounting only for the most physi- 
cally important terms. We will try to expose the approx- 
imations in a clear and logical way, providing the phys- 
ical justification as we go along. Resulting second-order 
model consist of nine coupled equations for the mean 
velocity and temperature gradients, four components of 
the Reynolds stresses, two components of the tempera- 
ture fluxes and the temperature variance. Thanks to the 
achieved simplicity of the model we found an approxi- 
mate analytical solution of these equations, expressing 
all nine correlations as functions of only one governing 
parameter, £(z)/L, where £{z) is the outer scale of tur- 
bulence (depending on the elevation z and also known as 
the "dissipation scale") and L - is the Monin-Obuckhov 
length. 

We would like also to stress, that in our approach 
£{z)/L is an external parameter of the problem. For 
small elevations z <C L, it is well accepted that £(z) is 
proportional to z, while the £(z) dependence is still under 
debate for z comparable or exceeding L. For z > L the 
assignment and discussion of the actual dependence of 
the outer scale of turbulence, £(z), which is manifested 
in the nature is out of the scope of this paper, and is 
remained for future work. At time being, we can analyze 
consequences of our approach for the following versions 
of £{z) dependence at z ^> L: 

• function £{z) is saturated at some level of the order of 
L. For concreteness we can take 

l/£(z) = ^(d lZ )- 2 + (d 2 L)- 2 , di~d2~l. (0.4) 

• £{z) is again proportional to z for elevations much 
larger than L: £{z) = d 3 z but with the proportionality 
constant d 3 < d\ . If so, we can also study the case £{z) » 
L even though such a condition may not be realizable in 
nature. In that case our analysis of the limit £{z) S> L 
has only a methodological character: it allows to derive 
an approximate analytic solution for all the objects of 
interest as functions of £(z)/L that, of course, is also 
valid for the outer scale of turbulence not exceeding a 
value of the order of L. 

It should be noticed that traditional turbulent closures 
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(including ours) cannot be applied for strongly strati- 
fied flows with Rig ra d ^ f (may be even at Ri gra d ~ 1)- 
The main reason is that these closures are roughly jus- 
tified for developed vortical turbulence, in which the 
eddy-turnover time is of the order of its life time; in 
other words, there are no well defined "quasi-particles" 
or waves. This is NOT the case for stable stratification 
with Rigrad f , in which the Brunt- Vaisala frequency 



N= ^f3dQ(z)/dz, (0.5) 

is larger then the eddy-turnover frequency 7. ft means 
that for Ri g r a d ^ 1 there are weakly decaying Kelvin- 
Helmholts internal gravity waves (with characteristic fre- 
quency N and decay time above I/7), propagating on 
large distances, essentially effecting on TBL, as pointed 
out by Zilitinkevich, (2002). In contrast to ambitious 
ZEKR attempt to describe entire TBL without limita- 
tion in values of Rigrad, and without explicit accounting 
for the internal gravity waves, we concentrate in our pa- 
per on self-consistent description of the lower part of the 
atmospheric TBL, in which turbulence has vortical char- 
acter and consequently, large values of Ri gra d do not ap- 
pear. We relate large values of Ri gra d in the upper part 
of TBL with contributions of the internal gravity waves 
to the energy and the energy flux in TBL, to the momen- 
tum flux, and to the production of (vortical) turbulent 
energy. Due to their instability in a shear flow, the waves 
can break and create turbulent kinetic energy. All these 
effects are beyond the framework of our paper. Their 
description in the upper "potential-wave" TBL and in- 
termediate region with the combined "vortical-potential" 
turbulent velocity field is in our nearest agenda. 

To make the paper more transparent for wide audience, 
not necessarily experts in atmospheric TBL, we attempt 
to present the material in a self-contained manner, and 
organized it as follows. 

In Sect. I A we used the Oberbeck-Boussinesq approx- 
imation and applied the standard Reynolds decomposi- 
tion (into mean values and turbulent zero-mean fluctu- 
ations of the velocity and temperature fields) to derive 
equations for the mean values and balance equations for 
all the relevant second-order correlation functions. In 
Sect. IB we demonstrate that the resulting balance equa- 
tions exactly preserve (in the non-dissipative limit) the 
total mechanical energy of the system, that consists of 
the kinetic energy of the mean flow, kinetic energy and 
potential energy of the turbulent subsystem. 

In Sect. II we describe the proposed closure procedure 
that results in a model of stably stratified TBL, that 
accounts explicitly for all the relevant second-order cor- 
relations. The third order correlations which appear in 
the theory are modeled in terms of second-order corre- 
lations in Sects. IIA and B. Further simplifications are 
presented in Sects. IIC and D for stationary turbulent 
flows in a plane geometry outside the viscous and buffer 
layers. In Sect. HE we suggest some generalization of 
the standard "wall-normalization" to obtain the model 



equations in dimensionless form with only one governing 
parameter, £(z)/L. 

Section III A contains approximate analytical solution 
of the model. It is shown that the analytical solution 
deviates from the numerical counterpart in less than a 
few percent in the entire interval < (i/L) < 00. 

The last Sect. Ill is devoted to a detailed description 
of our results: profiles of the mean velocity and potential 
temperature (Sect. IIIA), profiles of the turbulent kinetic 
and "temperature" energies, profiles of the anisotropy of 
partial kinetic energies (Sect. IIIB), profiles of the tur- 
bulent transport parameters u T and % T , profiles of the 
gradient- and flux-Richardson numbers Ri gra d and Rifl ux , 
and the dependence of the turbulent Prandl number Pr T 
vs. i/L and Ri gra d, Sect. IIIC. In conclusion Sect. HID, 
we consider the validity of the down-gradient transport 
concept (0.3) and explain why it is violated in the up- 
per part of TBL. The problem of critical Ri gra d is also 
discussed. 

In Appendix A we revisit the derivation of the 
Oberbeck-Boussinesq approximation. The reason for do- 
ing so stems from the fact that usual derivations are too 
specialized either to fluids whose density variation is very 
small, or to ideal gases. Thus it cannot be directly ap- 
plied, for example, to humid air. 

The detailed procedure of obtaining an approximate 
analytical solutions for all the objects of interest as func- 
tions of i/L is presented in Appendix B. Exact solu- 
tions are obtained for two cases: neutral stratification, 
i/L = 0, and (methodological) limit i/L ~ > 00. Both so- 
lutions arc corrected up to the first order in correspond- 
ing small parameters. The desired approximate analyti- 
cal solutions interpolate these two limits over all values 
oil/L. 

Appendix C contains a discussion of the applicability 
for the stratified flows of applied closure of the triple 
correlations via second order ones with the help of di- 
mensional reasoning. 

Appendixes A-C are available in the online version of 
this paper or at Los- Alamos archive # nlin.CD/06f 201 8 



I. SIMPLIFIED DYNAMICS IN A STABLY 
TEMPERATURE-STRATIFIED TBL AND THEIR 
CONSERVATION LAWS 

The aim of this section is to consider the simplified 
dynamics of a stably temperature-stratified turbulent 
boundary layer, aiming finally at an explicit description 
of the height dependence of important quantities like the 
mean velocity, mean temperature, turbulent kinetic and 
potential energies, etc. In general one expects very differ- 
ent profiles from those known in standard (unstratified) 
wall-bounded turbulence. We want to focus on these dif- 
ferences and propose that they occur already relatively 
close to the ground allowing us to neglect (to the leading 
order) the dependence of the density on height and the 
Coriolis force. We thus begin by simplifying the hydro- 
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dynamic equations which are used in this section. 



A. Simplified hydrodynamic equations and 
Reynolds decomposition 

First we briefly overview the derivation of the govern- 
ing equations in the Boussinesq approximation. For the 
extensive details the reader is referred to the Appendix 
A. 

The system (Al) of hydrodynamic equations describ- 
ing a fluid in which the temperature is not uniform con- 
sists of the Navier-Stokes equations for the fluid velocity, 
U(r,t), a continuity equation for the space and time de- 
pendent (total) density of the fluid, p(r,t), and of the 
heat balance equation for the (total) entropy per unit 
mass, S(r,t), Landau and Lifshitz, 1987. 

These equations are considered with boundary condi- 
tions that maintain the solution far from the equilibrium 
state, where TA — S = 0. These boundary conditions are 
U = at zero elevation, U = const at a high elevation of 
a few kilometers. This reflects the existence of a wind at 
high elevation, but we do not attempt to model the phys- 
ical origin of this wind in any detail. The only important 
condition with regards to this wind is that it maintains a 
momentum flux towards the ground that is prescribed as 
a function of the elevation. Similarly, we assume that a 
stable temperature stratification is maintained such that 
the heat flux towards the ground is prescribed as well. 

We neglect the viscous entropy production term as- 
suming that the temperature gradients are large enough 
such that the thermal entropy production term domi- 
nates. For simplicity of the presentation we restrict our- 
selves by relatively small elevations and disregard the 
Coriolis force (for more details, see Wyngaard, 1992). 
On the other hand we assume that the temperature and 
density gradients in the entire turbulent boundary layer 
are sufficiently small to allow employment of local ther- 
modynamic equilibrium. In other words, we assume the 
validity of the equation of state. 

As a "basic reference state" (BRS) denoted hereafter 
by a subscript 'V we use the isentropic model of the 
atmosphere, where the entropy is considered space ho- 
mogeneous. Now assuming a smallness of deviations of 
the density and pressure from their BRS values and ex- 
ploiting the equation of state, one obtains a simplified 
equation (A7) from the complete system (Al), which is 
already very close to the standard Navier-Stoks equation 
in the Boussinesq approximation. Introducing (general- 
ized, see eq. (A9)) potential temperature, one results in 
the well-known system (A12) of hydrodynamic equations 
in the Boussinesq approximation. 

We confine ourselves to regions which are close to 
the ground hence neglect the dependence of the density 
on height, thus we replace V(p/pb) => (^p)/Pb, and 
V(pbU) => PbV • U. The resulting equations read 



vu v P _ A7/ ve d 

Vt pb T>t 



= X A9 d . (1.1) 



Here V/CDt = d/dt+U- V is the convection time deriva- 
tive, p - pressure, pb is the density in BRS, /3 = g(3 is 
the buoyancy parameter (g is the gravity acceleration 
and (3 is the thermal expansion coefficient Eq. (A6d), 
which is equal to 1/T, reciprocal molecular tempera- 
ture, for an ideal gas), 6d is the deviation of the po- 
tential temperature from its BRS value, v - dynamical 
viscosity and \ is the- dynamical thermal conductivity. 
To develop equations for mean quantities and correla- 
tion functions one applies the Reynolds decomposition 

u = u + u, (u) = u, (u) = o,e d = e + e, (e d > = 

6 , (0) = , p = (p) + p , (p) = 0. Here the aver- 
age (• • • ) stands for averaging over a horizontal plane at 
a constant elevation. This leaves the average quantities 
with a z,t dependence only. Substituting in Eqs. (1.1) 
one gets equations of motion for the mean velocity and 
mean temperature profiles 



DUi 
Dt 



+ V,-Ti. 



VM -A6, ^ + V.F = 0. (1.2) 



Pb 



Dt 



Here D/Dt = d/dt+U • V is the mean convection deriva- 
tive. The total (molecular and turbulent) momentum 
and thermal fluxes are 



-1/Vj Ui + T t j 



-xve + F, 



(1.3) 



where 7v. 



is the Reynolds stress tensor describ- 



ing the turbulent momentum flux, and F = (u9) is the 
turbulent thermal flux. In order to derive equations for 
these correlation functions, one considers the equations 
of motion for the fluctuating velocity and temperature: 



Du/Dt = — ii ■ VU - u ■ Vu 

-(Vp/p b ) + vAu-f3e, 
D6/Dt = -u- VQ-u- V6 + X A9 



u-Vu) (1.4a) 
0, 

u • V6) . (l-4b) 



The whole set of the second order correlation functions 
includes the Reynolds stress, r^ , the turbulent thermal 
flux, F, and the "temperature energy" Eg = (Q 2 ) /2, 
which is denoted and named by analogy with the tur- 
bulent kinetic energy density (per unit mass and unit 
volume), E K = \ (|w| 2 )/2 = Tr{r ij }/2. Using (1.4) one 
gets the following "balance equations" : 



Dnj d 



DFi 

UT + 

DE e 
Dt 



_d_ 

dXj 

+ e + V-T = 



Cij H~ T^ij ■> (1.5a) 
Ai + Bi, (1.5b) 
-F-V0. (1.5c) 



Here we denoted the dissipations of the Reynolds-stress, 
heat-flux and the temperature energy by 



£ ij — 



2 / duj duj 
\dx k dx k 

X<|V0| 2 ), 



e ■ = (v + )( — — 

1 ~ \ dx k dx k 



(1.6) 
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The last term on the LHS of each of Eqs. (1.5) describes 
spatial flux of the corresponding quantity. In models of 
wall bounded unstratified turbulence it is known that 
these terms are very small almost everywhere. We do not 
have sufficient experience with the stratified counterpart 
to be able to assert that the same is true here. Neverthe- 
less, for simplicity we are going to neglect these terms. 
It is possible to show that the accounting for these terms 
does not influence much the results. Note that keep- 
ing these terms turns the model into a set of differential 
equations which are very cumbersome to analyze. This 
is a serious uncontrolled step in our development, so we 
cross our fingers and proceed with caution. Since these 
terms are neglected we do not provide here the explicit 
expressions for Tijk, Ty, and T. 

The first term on the RHS of the balance Eq. (1.5a) 
for the Reynolds stresses is the Energy Production ten- 
sor Vij , describing the production of the turbulent kinetic 
energy from the kinetic energy of the mean flow, propor- 
tional to the gradient of the mean velocity: 

Vij = -r lk dUj/dxk - Tj k dUi/dxk ■ (1.7a) 

The second term on the RHS of Eq. (1.5a), dj, will be 
referred hereafter to as the "Energy conversion tensor" . 
It describes the conversion of the kinetic turbulent energy 
into potential energy. This term is proportional to the 
buoyancy parameter (3 and the turbulent thermal flux F: 

dj^-^Fidjz + FjSiz) . (1.7b) 

The next term in the RHS of Eq. (1.5a) is known as the 
"Pressure-rate-of-strain" tensor: 

Tlij = (psy/pb) , Sij = dui/dxj + duj/dxi . (1.7c) 

In incompressible turbulence its trace vanishes, therefore 
IZij does not contribute to the balance of the kinetic en- 
ergy. As we will show in Sec. II A, this tensor can be pre- 
sented as the sum of three contributions (Zcman, 1981), 

1l i j = R?J + R% + R%, (1.7d) 

in which R™ is responsible for the nonlinear process of 
isotropization of turbulence and is traditionally called the 
"Return-to-Isotropy" , R-J is similar to the energy pro- 
duction tensor (1.7a) and is called "Isotropization of pro- 
duction". New term, appearing in the stratified flow R-j, 
is similar to the energy conversion tensor (1.7b) and will 
be refereed to as the "Isotropization of conversion" . 

Consider the balance of the turbulent thermal flux F, 
Eq. (1.5b). The first term in the RHS, A describes the 
source of F and, by analogy with the energy-production 
tensor, Vij, is called "Thermal-flux production vector". 

Like Vij, Eq. (1.7b), it has the contribution, A*" , pro- 
portional to the mean velocity gradient: 

su se ,E0 su 

Ai = At +A t +A t , A % =-F-VUi, 
AT = -nj dQ/dxj , A™ = 2{iE e 5 lz , (1.7e) 



and two additional contributions, related to the temper- 
ature gradient and to the "temperature energy", Eg, and 
the buoyancy parameter. One sees, that in contrary to 
the oversimplified assumption (0.3b) the thermal flux in 
turbulent media cannot be considered as proportional to 
the temperature gradient. It has also a contribution, pro- 
portional to the velocity gradient and even to the square 
of the temperature fluctuations. Moreover, the RHS of 
the flux-balance Eq. (1.5b) has an additional term, the 
"Pressure-temperature-gradient vector" which, similarly 
to the pressure-rate-of-strain tensor (1.7d), can be di- 
vided into three parts (Zeman, 1981): 

B= (pVe/p b ) = B T<n + B sv + B Ee . (1.7f) 

As we will show in Sec. II A the first contribution, Bf° oc 
(u u Vj 0) is responsible for the nonlinear flux of F in the 
space of scales, toward smaller scales, similarly to the 
correlation (uuu), which is responsible for the flux of 
kinetic energy (u 2 ) toward smaller scales. The correla- 
tion Bf D oc (uuVi 9) may be understood as the nonlin- 
ear contribution to the dissipation of the thermal flux. 
Correspondingly we will call it "Rcnormalization of the 
thermal-flux Dissipation" and will supply it with a su- 
perscript " RD " . The next two terms in the decomposi- 
tion (1.7f) arc Bt oc S u and Bf oc Eg. They describe 
the renormalization of the thermal-flux production terms 
At oc S v and oc Eg, accordingly. 

B. Conservation of total mechanical energy in the 
exact balance equations 

The total mechanical energy of temperature stratified 
turbulent flows consists of three parts with densities (per 
unit mass): E = Ex+E K + E P , where E K = \U\ 2 /2 is the 
density of kinetic energy of the mean flow, E K = Ta/2 is 
the density of turbulent kinetic energy and E P = [3Eq/Sq 
is the density of potential energy, associated with turbu- 
lent density fluctuation p = f39pb, caused by the (poten- 
tial) temperature fluctuations 6, and S@ = d&/dz. Note 
that the formula for E-p appears different from the plane 
average of the £ P in Eq. (A14). In Appendix A 4 we 
show that in fact the difference between the two objects 
is a time independent total potential energy in the basics 
reference state, and therefore it can be considered as the 
zero level of the potential energy. 

The balance Eq. for Ex follows from Eq. (1.2): 

DE K /Dt+v(VjUif+Vj {U t n 3 ) = [source Erf+ryVj U t , 

(1.8a) 

with the help of identity: UiS/j Tij = Vj (UiTij) — TijVj Ui 
and definition (1.3). The terms on the LHS of this Eq., 
proportional to v and fij respectively, describe the dissi- 
pation and the spatial flux of Ex.- The term [source Ex] 
on the RHS of Eq. (1.8a) describes the external source 
of energy, originating from the boundary conditions de- 
scribed above and t^X^ Ui, describes the kinetic energy 
out-flux from the mean flow to turbulent subsystem. 
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The balance Eq. for the turbulent kinetic energy fol- 
lows directly from Eq. (1.5a): 

DE K /Dt+[e ii +V j T iij ]/2 = -T ij V j U i -PF z . (1.8b) ^ =Pb 

On the LHS of Eq. (1.8b) one sees the dissipation and 
spatial flux terms. The first term on the RHS originates 
from the energy production, \Vu, defined by Eq. (1.7a). 
This term has an opposite sign to the last term on the 
RHS of Eq. (1.8a) and describes the production of the 
turbulent kinetic energy from the kinetic energy of the 
mean flow. The last term on the RHS of Eq. (1.8b) origi- 
nates from the energy conversion tensor | Cu, Eq. (1.7b), 
and describes the conversion of the turbulent kinetic en- 
ergy into potential one. 

According to the last of Eqs. (1.5), one gets the bal- 
ance equation for the potential energy E P : multiplying 
Eq. (1.5c) for Eg by (3/S&: 



DE P /Dt + f3 



/S e = I3F Z 



(1.8c) 



The RHS of this Eq. [coinciding up to a sign with the last 
term on the RHS of Eq. (1.8b)] is the source of potential 
energy (from the kinetic one). 

In the sum of the three balance equations the conver- 
sion terms (of the kinetic energy from the mean to turbu- 
lent flows and of the turbulent kinetic to potential one) 
cancel and one gets the total mechanical energy balance: 

DE /Dt + [diss E] + V [AuxE] = [source E K \ . (1.9) 

This equation exactly respects the conservation of total 
mechanical energy in the dissipation- less limit, irrespec- 
tive of the closure approximations. This is because the 
energy production and conversion terms are exact and do 
not require any closures, while the pressure-rate-of-strain 
tensor, that requires some closure, does not contribute to 
the total energy balance. 



II. THE CLOSURE PROCEDURE AND THE 
RESULTING MODEL 

In this section we describe the proposed closure proce- 
dure that results in a model of stably stratified TBL. In 
developing this model we rely strongly on the analogous 
well developed modeling of standard (unstratified) TBL. 
The final justification of this approach can be only done 
in comparison to data from experiments and simulations. 
We will do below what we can to use the existing data, 
but we propose at this point that much more experimen- 
tal and simulational work is necessary to solidify all the 
steps taken in this section. 



tuating part of the pressure 
son's equation for p follows 



from 



The 
Eq. 



Pois- 
(1.4): 



j (uiUj - (uiUj) + UiUj + UjUi)+PV z 6 . 

The solution of this equation includes a harmonic part, 
Ap = 0, which is responsible for sound propagation 
and docs not contribute to turbulent dynamics at 
small Mach numbers. Thus this contribution can be 
neglected, the inhomogeneous solution includes three 
parts p= pb\Puu +PUu +Pe], where 



p uu = A *V 'jV j ((UiUj) - UiUj) , (2.1) 

PUu = A^ViVj {UiUj + UjUi) , p e = /3A" 1 V 2 9 , 



and the inverse Laplace operator A -1 is defined as usual 
in terms of an integral over the Green's function. 

Correspondingly the correlations TZij and 3 consist of 
three terms, Eqs. (l-7d) and (l-7f), in which 



pip _ 
R tj = 



(PUu ■ 



pic 



{pesij) ,(2.2) 



(p uu ve) , b = (puuve) , B i = (p e ve) . 



All of those terms originating from p uu are the most prob- 
lematic because they introduce coupling to triple corre- 
lation functions: Rfj cx {uiUjUk) and B RD cx (u 2 V#). 
Thus they require closure procedures whose justification 
can be only tested a-posteriori against the data. 

Having in mind to simplify the model in most possible 
manner, we adopt for the diagonal part of the Return- 
to-Isotropy tensor, the simplest Rota form (Rotta, 1951) 



Rl. 



-7w (r ri -2£ K /3) , 



(2.3a) 



in which 7 RI is the relaxation frequency of diagonal com- 
ponents of the Reynolds-stress tensor toward its isotropic 
form, 2E K /3. The parametrization of 7 M will be dis- 
cussed later. The tensor R^ 1 is traceless, therefore the 
frequency 7 M must be the same for all the diagonal com- 
ponents of Rfl ■ On the other hand there are no reasons 
to assume that off-diagonal terms have the same relax- 
ation frequency. Therefore, following L'vov et al. (2006) 
we assume that 



DM 

Rij 



-JmTi 



I] 7 



(2.3b) 



with, generally speaking, 7 RI ^ 7 M . Moreover, on intu- 
itive level, we can expect, that off-diagonal terms should 
decay faster then the diagonal ones, i.e. 7 RI > 7 RI . In- 
deed, our analysis of DNS results, see Eq. (B8) shows 
that 7ri/7 ri ~ 1.46. 

The term B UD also describes return-to-isotropy due to 
nonlinear turbulence self interactions (Zcman, 1981), and 
may be modeled as: 



\. Pressure-Rate-of-Strain tensor TZij and 
Pressure-Temperature-Gradient vector B 



The correlation functions TZij and 23, de- 
fined by Eqs. (1.7c) and (l-7f), include flue- 



Br = -j RD Fi 



(2.3c) 



This equation dictates the vectorial structure of Bf D cx 
Fi, which will be confirmed below. The rest can be un- 
derstood as the definition of the 7rd as the relaxation 
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frequency of the thermal flux. Its parametrization is the 
subject of further discussion in Sec. II D. 

The traceless "Isotropization-of-Production" tensor, 
R-J, has a very similar structure to the production ten- 
sor, Vij, Eq. (1.7), and thus is traditionally modeled in 
terms of Vij (Pope, 2001): 

RIJ ~ -C IP (Vij - % V/3) , V = Tt {Vij} . (2.3d) 

The accepted value of the numerical constant C IP = | 
(Pope, 2001). 

The traceless "Isotropization-of- Conversion" tensor, 
RJf does not exist in unstratificd TBL. Its structure 
is very similar to the conversion tensor, Cy, Eq. (1.7). 
Therefore it is reasonable to model it in the same way in 
terms of (Zeman, 1981): 

R% ~ -C 1C {dj - % C/3) , C = Tr {C v } , (2.3e) 

with some new constant C IC . 

The renormalization of production terms B i and B i 
are very similar to the corresponding thermal flux pro- 
duction terms, A i and A i , defined by Eqs. (1.7). 
Therefore, in the spirit of Eqs. (2.3d) and (2.3e), they 
are modelled as follows: 

B- U = (C su - l)Af - (1 - C SU )(F -V)Ui, (2.3f) 

E8 EQ 

Bi =-(C Ee +l)A =-20(C B ,+l)E o 5 ix .(2.3g) 
Using this and (2.1) one finds the sign of C Ee : 

-p{C Ee + l)E e = (p e V z 6)=f3((V z 6)A- 1 (V z 6)), 

c Ee = (i + ((v z 6)A-\v z e))/(e 2 )) < o . (2.4) 

To estimate C Ef) we assume that on the gradient 
scales the temperature fluctuations are roughly isotropic 
(ZEKR-paper), and therefore we can estimate A = + 
X7 y + ~ 3V^. Introducing this estimate and integrat- 
ing by parts leads to C E0 w —2/3. 

B. Reynolds-stress-, thermal-flux-, and 
thermal-dissipation 

Far away from the wall and for large Reynolds num- 
bers the dissipation tensors are dominated by the viscous 
scale motions, at which turbulence can be considered as 
isotropic. Therefore, the vector e should vanish, while 
the tensor Sij, Eq. (1.7), should be diagonal: 

e» = 0, Sij =2^ uu E K 5ij/3, (2.5a) 

where the numerical prefactor | is chosen such that j uu 
becomes the relaxation frequency of the turbulent kinetic 
energy. Under stationary conditions the rate of turbu- 
lent kinetic energy dissipation is equal to the energy flux 
through scales, that can be estimated as (uuu) /£, where 
I is the outer scale of turbulence. Therefore, the natural 
estimate of j uu involves the tripple-velocity correlator, 



luu ~ ({uuu) 1 1 (uu)), exactly in the same manner, as the 
Return-to-Isotropy frequencies, 7 RI and j K1 in Eqs. (2.3a) 
and (2.3b). Similarly, 

e = le eE e , m ~ (60u) / £ (60) . (2.5b) 

C. Stationary balance equations in plain geometry 

In the plane geometry, the equations simplify fur- 
ther. The mean velocity is oriented in the (stream- 
wise) x direction and all mean values depend on the 
vertical (wall-normal) coordinate z only: U = U(z)x, 
9 = 6», nj = Tij(z), F = F(z), E e = E e (z). 
Therefore (U ■ V) (...) = 0, and in the stationary case, 
when d jdt = 0, the mean convective derivative van- 
ishes: D /Dt = 0. Moreover due to the y — > — y sym- 
metry of the problem the following correlations vanish: 
T xy = T yz = F y = 0. The only non-zero components of 
the mean velocity and temperature gradients are: 

S„ = dU/dz , S e = dO/dz . (2.6) 

Finally, the acceleration due to gravity g is directed ver- 
tically and thus the buoyancy parameter has just one 
component: (3 = z(3. 

1. Equations for the mean velocity and temperature profiles 

Having in mind Eqs. of Sec. II C and integrating Eqs. 
(1.2) for U x and 9 over z, one gets equations for the 
total (turbulent and molecular) mechanical-momentum 
flux, t(z), and thermal flux, F, toward the wall 

t xz (z) = -vS v + t xz => t xz (0) = -r* , (2.7a) 
F z (z) - - X S B +F Z ^ F z (0) = —F* . (2.7b) 

The total flux of the x-component of the mechanical mo- 
ment in z-dircction is pwf xz {z) = J dz(d (p) / dx) + const. 
Generally speaking, t xz (z) depends on z. For example, 
for the pressure driven planar channel flow (of the half- 
wight S) Pb T xz (z) = (d ( P ) /dx){6 - z) < 0. 

Relatively close to the ground, where z <C 5, the z de- 
pendence of t xz (z) can be neglected. In the absence of 
the mean horizontal pressure drop and spatial distributed 
heat sources r and F are z-independent, and thus equal 
to their values at zero elevation, as indicated in Eqs. (2.7) 
after "=>"-sign. Notice, that in our case of stable strati- 
fication both vertical fluxes, the x-component of the me- 
chanical momentum, t XZi and the thermal flux, F z , are 
directed toward the ground, i.e. negative. For the sake 
of convenience, we introduce in Eqs. (2.7) notations for 
their (positive) zero level modulus: t* and F*. 

Recall that in the plain geometry U z = 0. Nevertheless 
one can write the equation for U z : 

d(T zz + (p)/ Pb )/dz = l3e, (2.8) 
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which describes a turbulent correction (oc t zz ) to the hy- 
drostatic equilibrium (A4b). Actually, this equation de- 
termines the profile of (p), that does not appear in the 
system of balance equations (2.7). 

2. Equations for the pair (cross)-correlation functions 

Consider first the balance Eqs. (1.5a) for the diago- 
nal components of the Reynolds-stress tensor in algebraic 
model (which arises when we neglect the spatial fluxes) : 

TE K + 3 7ai T ra /2 = - (3 - 2 C IP ) t xz S v - C ia (3 F z , 
TE K + 3 7ri t to /2 = -C 1p t xz S v - Cio P F z , (2.9) 
TE K + 3j ri t zz /2 = -C ip t xz 

where T = 7tm — 7ru . The LHS of these equations in- 
clude the dissipation and Return-to-isotropy terms. On 
the RHS we have the kinetic energy production and 
isotropization of production terms (both proportional to 
S v ) together with the conversion and isotropization of 
conversion terms, that are proportional to the vertical 
thermal flux F z . The horizontal component of the ther- 
mal flux, F x , does not appear in these equations. 

Equations (2.9) allow to find anisotropy of the 
turbulent-velocity fluctuations and to get the balance 
Eqs. for the turbulent kinetic energy with the energy 
production and conversion terms on the RHS: 

3r xx = 2{[2(l-C IP )r„ u / 7RI + l]£ K (2.10a) 
-{3-2C IP +C lc )0F z / lm }, 

3t to = 2{[(Gp-1)IW 7 r I + 1]£k (2.10b) 
-(C 1P + C IC )0F z / 7m }, 

Zt zz = 2{[(C IP -l)r uu /7 RI + l]£; K (2.10c) 
-(C ip -2C ic -3)/3F z / 7ri }, 
T UU E K = - Txz S v +(}F z , (2.10d) 

Equation (2.10d) includes the only non-vanishing tan- 
gential (off-diagonal) Reynolds stress t xz and has to be 
accompanied with an equation for this object: 

(C IP - 1)t zz S v + (1 + C lc )(3F x . (2.10e) 

This equation manifests that the tangential Reynolds 
stress t xz , that determines the energy production [ac- 
cording to Eq. (2.10d)], influences, in its turn, on the 
value of the streamwise thermal flux F x , which therefore 
effects on the turbulent kinetic energy production. 

As we mentioned, in the plain geometry F y = 0. Equa- 
tions (1.5b) for the F x and F z in this case take the form: 

lunF x = -(t xz S b +C sv F z S u ) , (2.11a) 
lnuF z = - (r zz S e + 2 C Ee 0E e ) , (2.11b) 

in which the RHS describes the thermal-flux production, 
corrected by the isotropization of production terms. 



The last Eq. (1.5c) for Eg, represents the balance be- 
tween the dissipation (LHS) and production (RHS): 

lee E e = -F z S e . (2.11c) 

D. Dimensional closure of time-scales and the 
balance equations in the turbulent region 

At this point we follow tradition in modeling of all the 
nonlinear inverse time-scales by dimensional estimates 
(Kolmogorov, 1941): 

Tim — c U u \/ E K J £ , 7ri = C Rl7mi , (2.12) 
7 ri = C ri7ri , 7 ee = Ceeluu , 7 rd = C u gj uu . 

Remember that t is the "outer scale of turbulence" . This 
scale equals to z for z < L, where L is the Monin- 
Obukhov length (definition is found below). 

Detailed analysis of experimental, DNS and LES data 
(see L'vov et al., 2006, and references therein) shows that 
for unstratified flows, g = 0, the anisotropic boundary 
layers exhibits values of the Reynolds stress tensor that 
can be well approximated by the values t xx — E K , T yy = 
t zz = E K /2. In our approach this dictates the choice 
C RI = 4(1 — Cip). Also we can expect that t vv is almost 
not affected by buoyancy. This gives simply C IC = — C IP . 
If so, Eqs. (2.10) with the parametrization (2.12) can be 
identically rewritten as follows: 

w ^ill- ^ Ek ■ & F ' 

T XX - & K , Tyy— , T ZZ — + , 

A fUU A Zj Zj j uu 

1uuE K = f3F z - t xz S u , j uu = c U u\[~E^ll, (2.13a) 
4 C M 7uu t X2 = (3 F x - t zz S v . 

For completeness we also repeated here the parametriza- 
tion (2.12) of 7uu . Finally we present the version of the 
balance Eqs. for the thermal flux (2.11a), (2.11b), and 
for the "temperature energy", (2.11c), after all the sim- 
plified assumptions: 

Cee luuEg — —F z S e , 

C u e luuF x = - ( Txz S e + C SV F Z S v ) , (2.13b) 
C U 9luuF z = - {t zz S q +2C Ee 0Ee) . 

E. Generalized wall normalization 

The analysis of the balance Eqs. (2.13) is drastically 
simplified if they are presented in dimcnsionlcss form. 
Traditionally, the conventional "wall units" are intro- 
duced via the wall friction velocity, w», and the viscous 
length-scale, A*, defined by u* = y/rZ, A* = v/u*. In ad- 
dition to w* , and A* , we introduce a new wall unit for the 
temperature, 0* = F^/u*, defined via the thermal flux at 
the wall and friction velocity. Subsequently, r + = r/A* , 

t + = t\*/u*, U + = U/ut, p + = p/pbul, 9 + = 9/0*, 
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t+ — F + 4- ft /9r 
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r yy k / ' 
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9 -r+ — P+ ^ 

^ ' ZZ — K 1 / C «W 
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li Li K 2 


^~XZ^u 5 




(2.17c) 






1 V *— ' r - « 


(2.17d) 


Cee c Mtt \J E^Eg = - 






(2.17c) 


C u8 c uu \[~E+F+ = - 


+ ct 


- C SU F+ S* , 


(2.17f) 


C u6 c U u\/e+F+ = - 


~ T zz s i - 


— 2 C Ee ^ £/g~ , 


(2.17g) 



6*+ = 0/0*, etc. Then the governing Eqs. (1.1) take the 
form: 

V+U + /Vt+ + V + p + = zQ+/L + + A+W+, 

V+e+/Vt+ = A+6+/Pr. (2.14) 

These Eqs. include two dimensionless parameters: the 
conventional Prandtl number Pr= v j 'k, and L + - the 
Monin-Obukhov length L measured in wall units: L = 
u l/PE r , L + = L/A*. Following ZEKR-paper, we used 
here the modern definition of the Monin-Obukhov length, 
which differs from the old definition by the absence of 
the von-Karman constant k w 0.436 in its denominator 
(Monin and Obukhov, 1954). 

Outside of the viscous sub-layer, where the kinematic 
viscosity and kinematic thermal conductivity can be ig- 
nored, L + is the only dimensionless parameter in the 
problem, which separates the region of weak stratifica- 
tion, z + < L+, and the region of strong stratification, 
where z + > L + . 

Given the generalized wall normalization we introduce 
objects with a superscript " + " in the usual manner: 

S+ = US V , S+ = A* SJ0. , 7 + = Ul , r+ = nj/ui , 
F+ = F/u*9* , E+ = Eg/91 . (2.15) 

In the turbulent region, governed by L + only, Eqs. (2.7) 
simplify to t+ z — — 1, F+ = — 1. Then, in the wall units 
Eqs. (2.13) can be reduced to Eqs. (Bl) for the six pro- 
files 5+, 5+, E+, Eg and F+ and r+, presented in the 
Appendix B 1. This equation can be effectively analyzed, 
see next Sec. Ill A. 



F. Rescaling symmetry and ^-representation 

Outside of the viscous region, where Eqs. (2.13) were 
derived, the problem has only one characteristic length, 
the Monin-Obukhov scale L. Correspondingly, one ex- 
pects that the only dimensionless parameter that governs 
the turbulent statistics in this region should be the ratio 
of the outer scale of turbulence, £, to the Monin-Obukhov 
length-scale L, which we denote as ft = £/L = £ + /L + . 
Indeed, introducing "J-objects" : 

& = t/L, sl = s^e+ , si = s+e+ , (2.16) 

and using Eqs. (2.15) one rewrites the balance Eqs. (2.13) 
as follows: 



These equations are the main result of current Sec. II. It 
may be considered as "Minimal Model" for stably strat- 
ified TBL, that respects the conservation of energy, de- 
scribes anisotropy of turbulence and all relevant fluxes 
explicitly and, nevertheless is still simple enough to al- 
low comprehensive analytical analysis, that results in ap- 
proximate analytical solution (with reasonable accuracy) 
for the mean velocity and temperature gradients S v and 
S e , and all second-order (cross)-correlation functions. 

As expected, the only parameter appearing in the Min- 
imal Model (2.17) is ft. The outer scale of turbulence, 
£, does not appear by itself, only via the definition of ft 
(B2). Therefore our goal now is to solve Eqs. (2.17) in or- 
der to find five functions of only one argument ft: S* , S^, 
E£, Eg and F+ . After that we can specify the depen- 
dence £ + {z + ) and then reconstruct the z + -dependence of 
these five objects. 



III. RESULTS AND DISCUSSION 

A. Analytical solution of the Minimal-Model 
balance equations (2.17) 

This Subsec. is devoted to analytical and numerical 
analysis of the Minimal-Model (2.17). Example of nu- 
merical solution of Eqs. (2.17) (with some reasonable 
choice of the phcnomcnological parameters) is shown in 
Fig. 1. Nevertheless, it would be much more instructive 
to have approximate analytical solutions for all correla- 
tions that will describe their ^-dependence with reason- 
able accuracy. The detailed procedure of finding these 
solutions is presented in Appendix B. A brief overview 
is as follows. 

Appendix B 1 presents the balance Eqs. (2.13) in gen- 
eralized wall units (2.15). Analysis of the resulting 
Eqs. (Bl) allows us to clarify the rescaling symmetry and 
to suggest " * " normalization (B2), in which the balance 
Eqs. (2.13) take final even simple form (2.17), which is 
the basis for further analysis. 

In Ap. B2 we show that Eqs. (2.17) can be reformu- 
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lated as polynomial Eq. (B6) of ninth order for the only 
unknown V E + . Analysis of its structure helps to formu- 
late effective interpolation formula (3.1), discussed below. 

Next, in Ap. B 3 we find the solutions (B7) of 
Eqs. (2.17) at neutral stratification, ft = 0, corrected 
by Eqs. (B12) up to linear order in ft. Its comparison 
with the existing DNS data results in an estimate for the 
constants C m ~ 1.46, and c uu 0.36, see L'vov at el. 
(2006) and references therein. 

Then, in Ap. B4b we considered the region ft — > oo. 
Even though such a condition may not be realizable in 
nature, from a methodological point of view, as we will 
see below, it enables to obtain the desired analytical ap- 
proximation. The ft ~ > oo asymptotic solution (B16) 
with corrections (B20), linear in the small parameter 
S = 2(c uu /ft) 4 / 3 is found. Now we are armed to sug- 
gest an interpolation formula 




that coincide with the exact solutions for ft = 0, 
Eqs. (B7), and for ft -> oo, Eq. (B16), including the 
leading corrections to both asymptotics, linear in ft, 

Eqs. (B12) and l/ft 4/3 , Eqs. (B20). Moreover, in the 
region ft <~ 1, Eq. (3.1a) accounts for the structure of ex- 
act polynomial Eq. (B6). As a result, the interpolation 
formula (3.1a) is close to the numerical solution with de- 
viations smaller than 3% in the entire region < ft < oo, 
see left middle panel on Fig. 1. Together with Eq. (B5b) 
it produces a solution for 5+, that can be written as 




where Lf = 3L+/14, = 3L + /11k and n is the 

von-Karman constant. This formula gives the same ac- 
curacy ~ 3%, see left upper panel in Fig. 1. We demon- 
strate below that the proposed interpolation formulae de- 
scribe the ^-dependence of the correlations with a very 
reasonable accuracy, about 10%, for any < ft < oo, see 
black dashed lines in Figs. 1. 

Unfortunately, a direct substitution of the interpola- 
tion formula (3.1) into the exact relation for S% obtained 
from the system (3.4) (see (B5c)) works well only for 
small ft, in spite of the fact that the interpolation for- 
mula is rather accurate in the whole region. The reason 
is a small denominator in Eq. (B5c) for large ft . We need 
therefore to derive an independent interpolation formula 
for S% . Using the expansions (B10) for small ft <C 1 and 
(B20) for large f ' > 1 we suggest 



si 



, + S+ +6( Cmi a) 4 / 3 S+~ 



(l + a£+/L+) i/3 



(3.1c) 



in which 

o+ _ 2 i/4 r ir 1 / i p+ 

St? = -C ue (2C m - (11 C u6 - 3 C su )/3 S+°°L + )/L+, 
St 00 = -14(C SU - 4C ue /3)/3L+ 
and a satisfies 

Stj + = Si x L + + 6Si?L + (c uu a)^ 3 -AaS e , /3, 
with 

= -°ue (3/4C RI - 22 + ZC sv /C ue ) /24C m . 

Equation (3.1c) is constructed such that the leading and 
sub-leading asymptotics for small and large ft coincide 
with first two terms in the exact expansions "almost" 
neural stratification, (B10), and extremely strong strati- 
fication, (B20). As a result Eq. (3.1c) approximates the 
exact solution with errors smaller then 5% for ft < 1 and 
ft > 50 and with errors smaller than 10% for any ft, see 
left upper panel in Figure 1. 

Substituting the approximate Eqs. (3.1) into the ex- 
act relations (B5d) and (B5e) and one gets approximate 
solutions Eq and F+ with errors smaller than 10%, see 
lower panels in Figure 1. 



B. Mean velocity and temperature profiles 

In principle, integrating the mean shear 5+ and the 
mean temperature gradient Si, one can find the mean 
velocity and temperature profiles. Unfortunately, to do 
so we need to know Si and Si as functions of the eleva- 
tion z , while in our approach they are found as functions 
of l/L. Remember, that the external parameter £ is the 
outer scale of turbulence that depends on the elevation 
z: I = £(z). For z<Lwe can safely take l—z, however 
for z > L the function £(z) is not found theoretically al- 
though it was discussed phenomenologically with support 
of observational, experimental and numerical data. It is 
traditionally believed that for z > L the scale £ saturates 
at some level of the order of L [see, e.g. Eq. (0.4)]. 

The resulting plots of U + are shown on Figure 2, up- 
per panel. Even taking £(z) — z one gets very similar 
velocity profile, see Figure 2, lower panel. With l(z) = z 
we find an analytical expression for the mean-velocity 
profile, using the interpolation Eq. (3.1b) for S+: 



U+(z) = - In 



zuo 1 + yi + (z/L 2 ) 



2/3 



Z 



(3.2) 
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FIG. 1: Color online. Log-log plots of the normalized velocity mean shears = £+ S+ and L+S+ (left upper panel) , normalized 

mean-temperature gradients S& = £ + and L + S& (right upper panel), the turbulent kinetic energy E+ and E+/ft 2/3 (left 

middle panel), partial kinetic energies Tu/E K (left middle panel), temperature energy E^ and E£ /ft 2 3 (left lower panel) and 

horizontal thermal flux F# and F£ /ft 2 ^ 3 (right lower panel) vs. ft — l/L = £ + / L + . Red and blue solid lines - exact numerical 
solutions before and after normalization by the large ft asymptotics, black dashed and dot-dashed lines - approximate analytical 
solutions. The region I > L may not be realized in the Nature. In this case it has only methodological character. 
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Here z„o is the roughness length. The ratios Lj/L are 
given with our choice of fitting constants. 

The resulting mean velocity profiles have logarithmic 
asymptotic for z < L and a linear behavior for z > L in 
agreement with meteorological observations. Usually the 
observations are parameterized by a so-called log-linear 
approximation: 



U + = k 1 ln(z/z u0 ) + z/Li . 



(3.3a) 



which is plotted in Figure 2 by dotted lines. One sees 
some deviation in the region of intermediate z. The rea- 
son is that the real profile [see, e.g. Eq. (3.2)] has a 
logarithmic term that saturates for z>L, while in the 
approximation (3.3a) this term continues to grow. To fix 
this one can use Eq. (3.2) (with L 2 — L\ for simplicity), 
or even its simplified version 



U+ = - In 



+ 



« z u0 VI + (z/Lrf U 



(3.3b) 



This approximation is plotted as a dashed line on Fig- 
ure 2 for comparison. One sees that our approxima- 
tion (3.3b) works much better than the traditional one. 
Thus we suggest Eq. (3.3b) for parameterizing meteoro- 
logical observations. 

The temperature profiles in our approach look similar 
to the velocity ones, see Figure 2, lower panel. They have 
logarithmic asymptotic for £ < L and linear behavior 
for £ > L. Correspondingly, they can be fitted by the 
log- linear approximation, like (3.3a), or even better, by 
improved version of it, like Eq. (3.3b). Clearly, the values 
of constant will be different: k => k t , Li => L\ t, etc. 



C. Profiles of second-order correlations 

The computed profiles of the turbulent kinetic and 
temperature energies, horizontal thermal flux profile and 
the anisotropy profiles are shown on Figure 1 in the mid- 
dle and lower panels. The anisotropy profiles, right mid- 
dle panel, saturate at £/L w 2, therefore they arc not 
sensitive to the z-dependence of £{z); even quantitatively 
one can think of these profiles as if they were plotted as 
a function of z/L. 

Another issue is the profiles of E+ (left middle panel) 
and of Eq and F+ (lower panels), that are oc (£/L) 2 / 3 for 
I 3> L (if available). With the interpolation formula (0.4) 
the profiles of the second order correlations have to sat- 
urate at levels corresponding to fi = 1. This sensitivity 
to the z-dependence of £{z) makes a comparison of the 
prediction with available data very desirable. 



D. Turbulent transport, Richardson and Prandtl 
numbers 



In our notations the turbulent viscosity and thermal 
conductivity, turbulent Prandtl number, the gradient- 



and flux-Richardson numbers are 

■*--£-s- ft(<,, £ 



Xt 

Pr T 

R-igrad 



F 1 , T+ 

l^ = ^ = C x (£t)^, 
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Pr T . 





(3.4a) 
(3.4b) 
(3.4c) 
(3.4d) 

(3.4e) 
(3.4f) 



With Eqs. (3.4a) and (3.4b) we introduce also two di- 
mensionless functions C v {fi) and C x {fi) that are taken 
as ^-independent constants in the down-gradient trans- 
port approximation (0.3) described in the Introduction. 
We will show, however, that the functions have a strong 
dependence on fi, going to zero in the limit fi — > oo 

as 1/fi 4 ^ 3 . Therefore this approximation is not valid for 
large fi even qualitatively. 



1. Approximation of down- gradient transport and its 
violation in stably stratified TBL 

As we mentioned in the Introduction, the concept of 
the down-gradient transport assumes that the momen- 
tum and thermal fluxes are proportional to the mean 
velocity and temperature gradients, see Eqs. (0.3): 



-v T S v 



F, 



"XtS'g ) 



(3.5) 



where v T and % T arc effective turbulent viscosity and 
thermal conductivity, that can be estimated by dimen- 
sional reasoning. Equations (0.3), giving this estimates, 
include additional physical arguments that vertical trans- 
port parameters should be estimated via vertical turbu- 
lent velocity, y 7 ?^, and characteristic vertical scale of 
turbulence, £ z . The relations between the scales £j in 
different j- directions in anisotropic turbulence can be 
found in the approximation of timc-isotropy, according 
to which 



7 



(3.6) 



Here 7 is a characteristic isotropic frequency of turbu- 
lence, that for concreteness can be taken as the kinetic en- 
ergy relaxation frequency 7„„. The approximation (3.6) 
is supported by experimental data, according to which in 
anisotropic turbulence the ratios £%/£j {i ^ j) are larger 



then the ratios L 



fru that are close to unity. 



With this approximations v T and % T can be estimated as 
follows: 



Z/ T = C V T ZZ /^ U 



Xt 



C x T zz /j u 



(3.7) 
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FIG. 2: Computed with Eq. (0.4) (for di — d,2 = 1) plots of U + (blue solid lines on left panels) and (red solid lines on right 
panels) vs \n(z/L) and vs. z/L (in inserts) for L + = 1000. In upper panels £(z) is taken from Eq. (0.4), while in lower panels 
take l{z) = z. Log-linear approximation (3.3a) to all profiles is shown by dotted lines, its improved version (3.3b) by dashed 
lines. The region I > L may not be realized in the Nature. In this case it has only methodological character. 



where, according to the approximation of down-gradient 
transport, the dimensionless parameters C„ and C x are 
taken as constants, independent of the level of stratifica- 
tion. 

In order to check how the approximation (3.5), (3.7) 
works in the stratified TBL for both fluxes, one can 
consider Eqs. (3.5) as the definitions of v T and Xt an d 
Eqs. (3.7) as the definitions of C v and C v . This gives 



7~xz 










r+ 9+ ' 

TzzD u 
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^Juu 




7~zz 
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' zz&e 



(3.8a) 
(3.8b) 



Recall, that in this paper the down-gradient approxima- 
tion is not used at all. Instead, we are using exact balance 
equations for all relevant second order correlation, includ- 
ing t xz and F x . Substituting our results in the RHS of 
the definitions (3.8) we can find, how C v and C x depend 
on V = £/L, that determines the level of stratification in 
our approach. 



The resulting plots of the ratios C v (ft)/C v (0) and 
C x (£*)/C x (0) are shown in the left-upper panel in Fig- 
ure 3. One sees that the C v (ft) and C x (fi) can be con- 
sidered approximately as constants only for £ < 0.2 L. 
For larger £/L both C v (£*) and C x (fi) rapidly decrease, 
more or less in the same manner, diminishing by an or- 
der of magnitude already for £ w 2 L. For larger £/L one 
can use the asymptotic solution (B13), (B16) for ^> 1 
according to which 

S+ ~ _ ~ ~ VjC/k ~ £ T ~ ft 2 ? 3 (c, q \ 

This means that both functions vanish as 

/r\4/3 /L\ i/3 

C v (i*) nO.Qllj) , C x (£*) ~ 0.003 ( -J ,(3.10) 

where numerical prefactors account for the accepted val- 
ues of the dimensionless fit parameters. 
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FIG. 3: Color online. Log-log plots of "down-gradient coefficient-functions" C„ (solid blue lines) and C x (red dashed lines) 
- upper panels; turbulent Prandtl number Pr T (green lines on middle panels) and Rifl ux (black solid lines), Ri gra d (black 
dashed lines) - on lower panels as function of ft = £/L (left panels) and vs. Ri gra d (right panels). Notice, that the presented 
dependencies have qualitative character, and the choice of constants C... depends on the actual functional form £(z). For 
simplicity, we took I (z) = z.The region £ > L may not be realized in the Nature. In this case it has only methodological 
character. 
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The physical reason for the strong dependence of C„ 
and C x on stratification is as follows: in the RHS of 
Eq. (2.10e) for the momentum flux and Eq. (2.11b) for 
the vertical heat flux there are two terms. The first ones, 
proportional to t zz and velocity (or temperature) gradi- 
ents correspond to the approximation (3.5), giving (in our 
notations) C v =const and C x =const, in agreement with 
the down-gradient transport concept. However, there are 
second contributions to the vertical momentum flux cx F x 
and to the vertical heat flux, that is proportional to f3Eg. 
In our approach both contributions are negative, giving 
rise the counter- gradient fluxes. What follows from our 
approach, is that these counter-gradient fluxes cancel (to 
the leading order) the down-gradient contributions in the 
limit ft — > oo. As a result, in this limit the effective tur- 
bulent diffusion and thermal conductivity vanish, mak- 
ing the down-gradient approximation for them (with con- 
stant C v and C x ) irrelevant even qualitatively for i > L. 

In our picture of stable temperature-stratified TBL, 
the turbulence exists at any elevations, where one can 
neglect the Coriolis force. Moreover, the turbulent ki- 
netic and temperature energies increase as (£/L) 2 / 3 for 
I > L, see left middle and lower panels in Figure 1. At 
the same time, the mean velocity and potential temper- 
ature change the (i?/L)-dependence from logarithmic lo 
linear, see Figure 2 and (modified) log-linear interpola- 
tion formula (3.3b). Correspondingly, the shear of the 
mean velocity and the mean temperature gradient sat- 
urate at some elevation (and at some i/L), and Ri gra d 
saturates as well. This predictions agree with large eddy 
simulation by Zilitinkevich and Esau (2006), see Fig. 5 
in ZEKR-paper, where Ri gra d can be considered as satu- 
rating around 0.4 for z/L w 100. 

Nevertheless, our analytical result of saturating Ri gra d 
disagree both with the ZEKR model and with various ob- 
servational, experimental and numerical data, collected 
in ZEKR-paper, see their Figs. 1, 2, upper panel of Figs. 
3 and 4, where various data are plotted as functions of 
the gradient Richardson number in the interval (0, 100). 
The conditions at which these data were obtained do not 
correspond to the the situation considered in this pa- 
per. Nevertheless, one cannot completely ignore the fact, 
shown in Fig.l of ZEKR-paper, that various data con- 
centrate around a linear dependence Pr T ~ Ri gra d in the 
two-decade interval 1 < Ri gra d < 100 (not withstanding 
the high degree of scatter). 

Notice that the turbulent closures of kind used above 
cannot be applied for strongly stratified flows with 
Ri g r a d ^ 1 (may be even at Ri gra d ~ 1). There are two 
reasons of that. The first one was mentioned in the In- 
troduction. Namely, for Ri gra d ^ 1 the Brunt- Vaisala 

frequency N = y/(3S s , N + = V l+' ' s l ar g er then the 



eddy-turnover frequency and therefore there are weakly 
decaying Kelvin-Helmholts internal gravity waves which, 
generally speaking has to be accounted in the momentum 
and energy balance in TBL. 

The second reason, that makes the results very sen- 



sitive to the contribution of internal waves follows from 
the fact that vortical turbulent fluxes vanish (at fixed 
velocity and temperature gradients). Therefore even rel- 
atively small contributions of a different nature to the 
momentum and thermal fluxes may be important. 

The final conclusion is that the TBL modeling at large 
level of stratification requires accounting for turbulence 
of the internal waves together with the vortical turbu- 
lence, and the analysis of available data calls for serious 
revision. Definitely, new observations, laboratory and 
numerical experiments with control of internal wave ac- 
tivity are very likely. 
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APPENDIX A: OBERECK-BOUSSINESQ 
APPROXIMATION AND CONSERVATION LAWS 

1. Basic hydrodynamic equations 

The system of hydrodynamic equations describing a 
fluid in which the temperature is not uniform consists of 
the Navier-Stokes equations for the fluid velocity, U(r, t), 
a continuity equation for the space and time depen- 
dent (total) density of the fluid, p(r,t), and of the heat 
balance equation for the (total) entropy per unit mass, 
S(r,t), Landau and Lifshitz, 1987 : 

VU 



vt 
dp 

dt 
VS 

} vt 

V 

vt 



-Vp-gp+V ■ pVU , 


(Ala) 


V-(pU) = 0, 


(Alb) 




(Ale) 




(Aid) 



Here V /Vt is the convective (Lagrangian) derivative, p is 
the pressure, g = —zg is the vertical acceleration due to 
gravity, \i and k are the (molecular) dynamical viscosity 
and heat conductivity. 

These equations are considered with boundary condi- 
tions that maintain the solution far from the equilibrium 
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state, where 14 = S = 0. These boundary conditions are 
U = at zero elevation, 1A — const at a high elevation of 
a few kilometers. This reflects the existence of a wind at 
high elevation, but we do not attempt to model the phys- 
ical origin of this wind in any detail. The only important 
condition with regards to this wind is that it maintains a 
momentum flux towards the ground that is prescribed as 
a function of the elevation. Similarly, we assume that a 
stable temperature stratification is maintained such that 
the heat flux towards the ground is prescribed as well. In 
the entropy balance Eq. (Ale) we have already neglected 
the viscous entropy production term, oc /i|VW| 2 , assum- 
ing that the temperature gradients are large enough such 
that the thermal entropy production term on the RHS of 
Eq. (Ale) dominates. Actually, this assumption is very 
realistic in our applications. For simplicity of the presen- 
tation we restrict ourselves by relatively small elevations 
and disregard in Eq. (Ala) the Coriolis force (for more 
details, see Wyngaard, 1992). 

On the other hand we assume that the temperature 
and density gradients in the entire turbulent boundary 
layer are sufficiently small to allow employment of local 
thermodynamic equilibrium. In other words, we assume 
the validity of the equation of state, and that the entropy 
S is a state function of the local values of the density and 
pressure: 

p = p(T,p), S = S(p,p). (A2) 

In the same manner we will neglect the temperature de- 
pendence of the dissipation parameters v and k. 

Pressure fluctuations caused by turbulent velocity fluc- 
tuations u propagate in a compressible medium with 
the sound velocity c s , causing time dependent density 
fluctuations of the order of (u/c s ) 2 po, where po is the 
mean density. Assuming that the square of the turbulent 
Mach number M 2 = (u/c s ) 2 is small compared to unity, 
we can neglect in Eq. (Alb) the partial time derivative: 
V-(joW) = 0, see e.g. Landau and Lifshitz, 1987. Even in 
tropical hurricanes of category five the mean wind veloc- 
ity U is below 300 Km/h. Usually, the turbulent velocity 
fluctuations u are less then C// 10, i.e. even in these ex- 
treme conditions u < 30 Km/h and M 2 < 10~ 3 (with 
c s ~ 1200 Km/h). Therefore the incompressibility ap- 
proximation V • {pU) = is well justified in atmospheric 
physics. In the ocean where the sound velocity is even 
larger and water velocities even smaller, this approxima- 
tion is quite excellent. 



2. Isentropic basic reference state 

In quite air, without turbulence, the pressure and the 
density depend on the elevation z simply due to grav- 
ity. For example, in full thermodynamic equilibrium the 
temperature is uniform, ^-independent, and the density 
decreases exponentially with the elevation. However, this 
equilibrium model of the atmosphere is not realistic, and 



cannot be used as a reference state about which the ac- 
tual dynamics is considered. A much better reference 
state is a situation in which the entropy is space homoge- 
neous. In this model the thermal conductivity (leading to 
the temperature homogeneity) is neglected with respect 
to heat transfer due to the vertical adiabatic mixing of 
air, leading to a z-independent entropy. Following tra- 
dition, we refer to the isentropic model as a "basic ref- 
erence state" and denote this state of the system with a 
subscript " t>": 

S b = S(p h ,p h ) = const, p b = p{T h ,p h ) . (A3) 

The first of Eq. (A3) relates the gradients of the pressure 
and density in this state: 

Another relation between pb and pb follows from the con- 
dition of hydrostatic equilibrium: 

Vp b = gpb ■ (A4b) 

Equations (A4) together with the first of Eqs. (A3) de- 
termine the density, pressure and temperature profiles in 
the isentropic basic reference state. 

3. Hydrodynamic equations in generalized 
Obereck— Boussinesq approximation 

a. Equations of motion 

Denote the deviations of the total density, pressure, 
temperature and entropy from the basic reference state 
as follows: 

P = p-Pb, p = p-pb, (A5) 
f = T-T h , S = S- S h ■ 

Following Obereck (1879) and Boussinesq (1903), assume 
that these deviations are small: p <C Pb, P <C Pb- 
Then one simplifies the full system of hydrodynamic 
Eqs. (Al) and rewrites them in terms of the fluid ve- 
locity U and the small deviations p and S (instead of p) . 
The first step is very simple: because of Eq. (A4b) 

Vp-gp=Vp-gp . (A6a) 

Next we should relate the deviations p, p and S. In the 
linear approximation Eq. (A2) yields: 

S = (dSb/dpb) p P+(dSb/dpb) p p ■ (A6b) 

With the help of Eqs. (A4) this gives: 

gp=^p-(3^T h S . (A6c) 
Pb c p 
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Here (3 = gf3is the buoyancy parameter, (3 is the thermal 
expansion coefficient and c p is the isobaric specific heat 
(heat capacity per unit mass at constant pressure) in the 
basic reference state: 



1 

Pb 



Pb T h 



dpb 
dS h 



(A6d) 



Now Eq. (A6a) in the linear approximation yields: 

'Pb 



Vp-gp = p h 



c 



yPb / >-p 

Then Eqs. (Ala) can be approximated as: 



(A6e) 



VU 

~vt 



-VI*. 

Pb 



f3 — S+— V-nVU. (A7) 

Cp Pb 



b. Generalized potential temperature 

To proceed, we generalize the notion of potential tem- 
perature which is traditionally defined as the temper- 
ature that a volume of dry air at a pressure p(z) and 
temperature T(z) would attain when adiabatically com- 
pressed to the pressure that exists at zero elevation 
z = 0. This potential temperature can be explicitly com- 
puted for an ideal gas with the result 



Q(z) = n( P ,/p(z)) 



( 7 -l)/7 



(A8) 



where 7 is the ratio of isobaric to isochoric specific heats, 
7 = c p /c v , and is the temperature at zero elevation. 

We want to generalize the notion of potential tempera- 
ture for an arbitrary stratified fluid requiring that in the 
isenotropic basic reference state it would be a constant 
6* = T*. A second requirement is that the definition 
will agree with Eq. (A8) for an ideal gas. Accordingly 
we define 



Q(z) = T* exp [(S(z) - 5b) /c p ] 



(A9) 



For more details see also Hauf and Holler (1987). Indeed, 
if we employ the equation of state and the equation for 
the entropy of an ideal gas, i.e. 



p = pT, 5 = In 



const 



(A10) 



one can easily check that Eq. (A8) is recaptured. 

For small deviations of 6 from the basic reference state 
value T„, Eq. (A9) gives up to linear order: 



(All) 



Now we can present Eqs. (Ala), (Alb) (with dpt/dt = 0, 
as explained) and (Ale) as follows: 

— = -V (j^J -(3Q d + isAU, (A12a) 



Vt 

V-(pbU)=Q 



Vt 



= xAe d . 



(A12b) 
(A12c) 



The dissipative terms are important only in the narrow 
region of the viscous sublayer, where we can safely ne- 
glect the z-dependence of pb, p and k, and consider the 
dynamical viscosity v = pjpb and dynamical thermal 
conductivity x = K /pb &s some z-independent constants. 

Note that for turbulence in liquids (water, etc.) one 
can simplify these equations further. There one can ne- 
glect the effect of adiabatic cooling (together with the 
compressibility), and simply use another reference state 
with constant temperature and density: 



T = T* , pb = p* , Pb = P* + gp*z . 



(A13) 



For this reference state standard reasoning (see, e.g. Lan- 
dau and Lifshitz, 1987) yields the same equations as 
Eqs. (A12) in which again (3, is given by Eq. (A6d) and 
it is a parameter characterizing a particular fluid. In 
this case pb = p*, independent of z, and the potential 
temperature = T, such that 6d — T d is a deviation 
of the total temperature T from its ground (bottom, or 
whatever) level T*. 

In our Eq. (A12) the situation is more general since 
we do not assume that reference state has the simple 
form (A13) (with pb = const). Importantly, on the 
RHS of Eq. (Al2a) the density Pb{z) is operated on 
by the gradient, and the buoyancy term — /39d involves 
9d 7^ Td, the deviation of the potential temperature de- 
fined by Eq. (All). This definition for liquids has noth- 
ing in common with the standard meteorological defini- 
tion (A8). Notice also that for an ideal gas (3 — 1/T, 
and Eq. (A12a) coincide with that suggested in the book 
Kurbatsky (2000) and used in ZEKR-paper. 

It is important to realize that the approximate Eqs. 
(A12) conserve exactly an approximate expression for the 
total mechanic energy of the system in the dissipation- 
less limit. Consider the sum of the kinetic, £ K , and 
the potential energy £ P (calculated in the basic reference 
state): £ t = £ K + £ P , 

£ K = J dr p b , £ P = J drp b (3 ■ r @ d . (A14) 

One can check by direct substitution that this sum of 
energies is conserved by of motion Eqs. (A12) when v = 

x = o. 



4. Potential energy in stratified TBL 

In this Appendix we show that the potential energy of 
a stratified turbulent flow £ P , Eq. (A14), can be presented 
as the sum of the (time-independent) potential energy of 
the basic reference state, £ P and a "turbulent" potential 
energy, associated with temperature fluctuations, £ P with 
the density per unit mass E P = PEq/Sq: 



I 



dr(p b (3E e /S e ) 



(A15) 
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Actually, we want to discuss this issue in a more general 
case, when the stratification is caused by some "internal" 
parameter of the fluid, £, not necessarily the potential 
temperature. It can be the salinity of water in a sea, 
the humidity of the air, the concentration of particles 
co-moving with the fluid as Lagrangian tracers, etc. 

In the general case then the equation for the potential 
energy of a stratified fluid has the form: 



= 9 J p(r 



)zdxdydz . 



(A16a) 



In the BRS the potential energy reaches its minimum 
value referred to as the basic potential energy, £ P : 



= 9 J Ph{z) 



z dx dy dz . 



(A16b) 



where p b (z) is the density profile in BRS. Clearly, in equi- 
librium the fluid density p b (z) decreases with the eleva- 
tion: dp b /dz < 0. 

In the turbulent state the density deviates from p b (z): 
p(r,t) = p b (z) + p(r,t) and the mean potential energy 
(£ P ) exceeds £ P . Now we can write 



(A16c) 



We compute the turbulent potential energy , £ P , in the 
case when the internal parameter £ (temperature, etc.) is 
co-moving with the fluid clement as a Lagrangian marker. 
In the Lagrangian approach we can consider p as the La- 
grangian marker and introduce the variable z(p, i), which 
is understood as an elevation of the fluid element with 
the density p. Noticing that zdz = ^dz 2 , and integrating 
Eq. (A16a) by parts with respect of z 2 , we can present 
(£ P ) in the Lagrangian approach as: 



(£p> = 



(KP>*)] 2 ) dxdydp 



(A17) 



As a result of turbulent motion, the elevation z(p,t) at 
given p fluctuates and can be decomposed into the mean 
and fluctuating parts: 

z(p,t) = z b {p)+z(p,t), (z{p,t))=Q. (A18) 



Substituting z 2 (p,t) 



2z b z + z 2 in Eq. (A17) we 



have three contributions to the potential energy The first 
one (originating from z 2 ) describes the basic potential 
energy, Eq. (A16b). The second contribution, which is 
linear in z, disappears because (z) = 0. The last one 
describes the turbulent potential energy: 

£» = -§/ (Ffo*)] 2 ) dxdydp . (A19) 

Relating the density fluctuations p, around the BRS den- 
sity profile p b (z) (caused by the fluctuations z) with z 



Sp^^lsz. 
dz 



(A20) 



and returning back to the Eulerian description in 
Eq. (A19) one has: 



dp b (z) 



dz 



dxdydz . 



(A21) 



Here we used the transformation formula, similar to 
Eq. (A20): dp = [dp (z) jd z] dz. Equation (A21) allows 
one to introduce a local density of turbulent potential 
energy per unit mass, 



E P 



such that 



/ 



dp b (z) 



dz 



-i{P 



p b E P dx dy dz . 



(A22a) 



(A22b) 



For the particular case of temperature stratification 

dp b (z) 



dz 



~dO 
dz 



P = Pb f36, (3 = g(3, 



(A23) 



where /? is the thermal expansion coefficient, related, as 
shown, to the buoyancy parameter (3. With Eqs. (A23) 
and definition of the temperature energy Eg = | ($ 2 )> 
Eq. (A22a) goes to E P — PEg/S®, as declared. 



APPENDIX B: ANALYSIS OF THE 
MINIMAL-MODEL BALANCE EQUATIONS 

1. Wall unites and ^-representation 

The Eqs. (2.13) in the wall units take the form: 

2 ' 



T+ 

' XX 


F+ 

- E+ z r+ 


T+ 

' ZZ 


K F+ 
2 + 2L+7+„' 


7+ E+ 


'xz'-'u ' 


4 Cri ltu T xz 


F+ 

' zz u u ' 


CgglLK 




Cue Juu F x 


= ~ T xz^t _ CsuFj~ 5+ 


CuOlLFt 


_ _+ q+ O C ™ E 



with 



F+ 



-1 



(Bla) 
(Bib) 

(Blc) 

(Bid) 
(Blc) 

(Blf) 
(Big) 



Outside of the viscous region, the problem has only 
one characteristic length, the Monin-Obukhov scale L. 
Correspondingly, one expects that the only dimensionlcss 
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parameter that governs the turbulent statistics in this 
region should be the ratio of the outer scale of turbulence, 
£, to the Monin-Obukhov length-scale L, which we denote 
as ft = £/L = £ + /L + . Indeed, introducing "J-objects" : 

fi=l/L, S t = S+£+ , St = S+£+ , (B2) 

one rewrites the balance Eqs. (Bl) as follows: 

r+ = E++ft/2c uu ^E+, r+=E+/2, (B3a) 

2 r+ =E+- ft/c uu ^E+, (B3b) 

CuuE+ 3/2 = ftF+ - t+S* , (B3c) 

4 C m c uu ^E+t+ = £tF+ - r+ St , (B3d) 

Cee c uu \/e+E+ = -F+ S* , (B3e) 

C u e cuu \[e+F+ = -t+S* - C SV F+ St , (B3f ) 

C U cuu \[~E+F+ = -t+SI - 2 C Ee £ X E+ , (B3g) 

2. Polynomial formulation of balance Eqs. (B3) 

In Sect. B4b we justify the simple choice 

C Ee ^-2C ee /3, (B4) 

which is used presently. With it the equations (2.17) 
[except of (2.17d)] yield expressions for S 1 * , S&, E# and 
F+ as functions of ft and X: 



X = ^jE+(£t), 

fit T?+ _ r y 3 



(B5a) 
(B5b) 



' XZ i 

suv = - " c2 € ueX !?L > (Bsc) 



= -- 

F+{£ X ) = 



3c uu X 3 + llftF+ ' 
CuuCeeX z 



(B5d) 



T+Sl + C sv St v F+ (b5c) 

Cum C u eX 



Substituting Eqs. (B5) into the remaining Eq. (2.17d) 
one gets an equation for E K which can be represented as 

the 9 th order polynomial for X = ^E+{ft): 

-12 c 2 uu Cue T+ 2 £tF+X 2 + (3 c uu X 3 + 11 ftF+) 

x {8 c 2 uu C m C u g t+ 2 X 2 - (c uu X 3 - £tF+) (B6) 

x [ (2 C su + Cue) £ % Ft + CuuCueX 3 ] } = . 

This equation can be solved numerically to find the func- 
tion X(ft). After that Eqs. (B5) allow to find all the 



required correlations as functions of ft , see red solid lines 
in Figure 1. 

The following subsections are dedicated to finding of 
an approximate analytical solutions for all correlations 
that will describe their ^-dependence with reasonable 
accuracy. First we will find in the next Sec. B 3 the so- 
lution of Eqs. (B5) and (B6) for ft = , corrected up 
to linear order in ft. Then, in Sec. B4b we will find the 
asymptotic solution for ft — > oo with corrections, linear 
in the small parameter 5 = 2(c uu /ft) 4 / 3 . 



3. Solution for neutral and weak stratification 

For neutral stratification, when ft = 0, Eq. (B6) triv- 
ially yields: 



X = X = (8C RI ) 1/4 , E+ = E+ fi = 2V / 2G^ . (B7a 
With Eqs. (B5) this gives: 



n Y3 



<7* - - 

J e,0 — 



2c uu C ue F+ 
X 



(B7b) 
(B7c) 
(B7d) 



CeeXo ' 

c,+ t xz^b,o + C su Sl F+ 

F xfi = „ • (B7e) 

Cuu ^uBM) 

Remembering that S\ = 1/k and using experimental 
values in the lag-law region E+ ss 3.42 and k w 0.436 
one finds: 



C B 



1.46, c u 



0.36 . 



(B8) 



Using the relationship between E* and k as suggested in 
L'vov et al. (2006), E^ = 18k 2 , one finds a relationship 
between these two constants, 



4\/2 CuuC'u 



(B9) 



which agrees with the known values (B8) with an accu- 
racy that is better then 1%. 

In the case of weak stratification one can develop a 
perturbative approach to Eq. (B6), using ft <C 1 as a 
small parameter, to develop a solution in the form of a 
Taylor series in powers of ft : 

X = X Q + X x ft + X 2 ft 2 + ... (BlOa) 
For our purpose only the linear term in ft is required: 



Xt 



-F: 



. 2 T+ z 2 C u e + C su Xq 
2 CuuCueXg 



(BlOb) 



where X is given by Eq. (B7a). Using now the exact 
relations (B5) with £?+ = X 2 one can immediately find 
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the linear terms in ft in the expansions for E+{ft), Sj r , 



Si, E+ and F+: 



E+ = E+ + E+ 1 ft + . 
sl = S* + S* /* + . . 

e> e,0 e,l 



where 

^K,l 
^1 

D e,l 
<1 



c £ +1/2 



^K,0 



2F+' 



3 c £ +5/2 



6 ^f^ _ 11 + 3-^" 



+ 2 

K,0 



^( 1 + 3 #2+^ 



V K,0 



2C„^ 



3 C S u 
2 C u e 

. 2 



(Blla) 
etc. (Bllb) 



(B12a) 
, (B12b) 
(B12c) 



C u e 



3 Csu 



ZE, 



+ 2 



E, 



+ 2 



-11 + ^ - ( B12d ) 



El 



(B12e) 



+ 



uu^uey -^k,o 
2c uu C ue ^E+~ 



4. Solution for strong gratification 

To analyze Eqs. (Bl) for large introduce new vari- 
ables for the second order correlations e K , eg, and / x , 
according to 



= (^) 2/3 , ^ = . - 



ee „t2/3 



In these variables Eqs. (Bl) take the form: 



>ft 



(B13) 



a. Asymptotic solution for ft — > oo 

In the limit — ► oo, one takes 5 — in Eqs. (B14) 
and gets an asymptotic solution (denoted by "°°") 

(B16a) 
(B16b) 
(B16c) 
(B16d) 
(B16e) 



^ = 


^uu \ 


-4—) , 

Cee J ' 


s+°° = 


^(' 


- 2— \ , 
Cee ) 


s+°° = 


-S+°°[ 


Csu + 2C U $ 


e? = 


L+S+ x 
CuuCee 


J X 


o Ceo 


L+S'+ oc 


Cee 





The physics of the problem requires that S^ S+, e K and 
eg be positively definite. This constraints the possible 
values of the fitting constants. With C ES < 0, Eq. (B4), 
we have: 

cuu > , Cee > , C u e > C su Cg e /2\C Ee | . (B17) 

Notice, that the maximal value of Rifl U x, [defined in 
Eq. (3.4c)], which we denote as Rif^max, corresponds 
to ft — > oo and equals to \/L + S^°°. Taking for 
Ceo = —2/3 according to the estimate C Eg w —2/3, 
and, following ZEKR-paper, accept Cee = 1, one gets 
R-if, max = 3/14 w 0.21 in agreement with the experimen- 
tal value (which is between 0.20 and 0.25). With the 
choice 

C E e = -2/3 , Cee = 1 • (B18) 
Equations (B16b) give the following predictions: 





= L+S+-1, 


(B14a) 


r p ) 2/3 

^UU C K ) 


= { C Uu e K — 1) ^'"'"S'u 


-2c uu / K , (B14b) 


Cee c uu ee 




(B14c) 


C u 9 C U u fx 


= L + (S~£ + C su S*u ) : 


(B14d) 


Cuu ^k) ^ 


= L + S^ {c uu e K — 1) 






+4C Ee c uu ee , 


(B14e) 



s+ c 



s+°° 



11 

3 Cuu 

14 



10.1, => E+ « 4.7 ^ 2/3 , (B19a) 



max Rifl 
4 



— wO.21, (B19b) 



3L+ ' 

) • (B19c) 
L+S+°°/c uu , (B19d) 
17.2, F+ w6.7^ 2/3 .(B19e) 



4L+5'+ oc 
3 c^ii 



with a small parameter S 



5 = 2^ 



Cuu \ 

it; 



(B15) 



appearing in the LHS of the equations for the vertical 
momentum and thermal fluxes. 



b. Expansion around ft — > oo 

The small parameter in Eqs. (B14) allows one to find 
the solution as Taylor series in 5 <C 1: 

S+ = S+°° + S^ 5 + St?5 2 + ... , etc. (B20a) 
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For our purposes we need to know only the linear terms 
in 5: 



C+oc _ C u 6 



(B20b) 
(B20c) 



s _ , rv, ( iic ue jc sv _ 2r . ^ (B2Qd) 



6,1 ~ L+ V 3S£°°L+ 

o+oo r + 

e e,i _ 



too _ si k k _ o u m 

Jx,i - C+o°r+ p ' 



(B20c) 
(B20f) 



APPENDIX C: ON THE CLOSURE PROBLEM 
OF TRIPLE CORRELATIONS VIA SECOND 
ORDER CORRELATIONS 

Let us look more carefully at the approximation (2.12), 
which is 



'Yuu — c uu \/ E K J t , ^f UJ — Cm^iuu J 

7ri = <?ri7m ) 76i6i = Ceeluu , 7rd 



(CI) 



The dimensional reasoning that leads to this approxima- 
tion is questionable for problems having a dimcnsionlcss 
parameter ft . Generally speaking, all "constants" c... and 
C... in Eq. (CI) can be any functions of ft. Presently we 



just hope that a possible ft dependence of these functions 
is relatively weak and does not affect the qualitative pic- 
ture of the phenomenon. 

Moreover, even the assumption (2.5a) that the dissipa- 
tion of the thermal flux is proportional to the thermal 
flux and the assumption (2.5b) that the dissipation of 
Eg, e oc Eg are also questionable. Formally speaking, one 
cannot guarantee that the tripple cross-correlator (9uu) + 
that estimates e + , can be (roughly speaking) decomposed 
like (u9) a/ (uu), i.e really proportional to F — (u9) as it 
stated in Eq. (2.5a). Theoretically, one cannot exclude 
the decomposition (Quu) <~ (uu) ^/ (06), i.e. a contri- 
bution to £ oc E K . Similarly, the dissipation e in the 
balance (1.5c) of Eg, that is determined by the correla- 
tor (2.5b), is oc (89u), as it follows from the decompo- 
sition (69u) ~ (99) yj(uu) and is stated in Eq. (2.5b). 
This correlator allows, for example, the decomposition 
(99u) ~ (9u) y/(9l), i.e. contribution to e oc F. This 
discussion demonstrates, that the situation with the dis- 
sipation rates is not so simple, as one may think and thus 
requires careful theoretical analysis that is in our agenda 
for future work. Our preliminary analysis of this problem 
shows that all fitting constants are indeed functions of fi . 
Fortunately, they vary within finite limits in the entire 
interval < fi < 00. Therefore we propose that the ap- 
proximations used in this paper preserve the qualitative 
picture of the phenomenon. Once again, the traditional 
down-gradient approximation does not work even qual- 
itatively because corresponding "constants" C v and C x 
vanish in the limit ft — > 00. 
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